% This practical will work with a single subject's data from an emotional
% faces experiment (Elekta Neuromag data). You can get the data from:
% https://sites.google.com/site/ohbaosl/practicals/practical-data/emotional-face-processing-elekta-neuromag-data
% 
% Work your way through the script cell by cell using the supplied dataset.
% As well as following the instructions below, make sure that you read all
% of the comments (indicated by %), as these explain what each step is
% doing. Note that you can run a cell (marked by %%) using the ?Cell? drop
% down menu on the Matlab GUI.    

%%%%%%%%%%%%%%%%%%
%% SETUP THE MATLAB PATHS
% make sure that fieldtrip and spm are not in your matlab path

global OSLDIR;
    
tilde='/home/mwoolrich';
osldir=[tilde '/Desktop/osl1.2.beta.12'];

tilde='/Users/woolrich';
osldir=[tilde '/homedir/matlab/osl1.2.beta.17'];

addpath(osldir);
osl_startup(osldir);

%%%%%%%%%%%%%%%%%%
%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS

testdir=[tilde '/Desktop'];
testdir=[tilde '/homedir/matlab/osl_testdata_dir'];

datadir=[testdir '/faces_subject1_data']; % directory where the data is
%datadir=[testdir '/faces_subject1_data_nosss']; % directory where the data is

workingdir=[datadir]; % this is the directory the SPM files will be stored in
%workingdir=[datadir]; % this is the directory the SPM files will be stored in

cmd = ['mkdir ' workingdir]; unix(cmd); % make dir to put the results in

% Set up the list of subjects and their structural scans for the analysis 
% Currently there is only 1 subject.
clear spm_files;

% set up a list of mris
structural_files{1}=[datadir '/structurals/struct1.nii'];

% set up a list of SPM MEEG object file names (we only have one here)
spm_files{1}=[workingdir '/spm8_meg1.mat'];
spm_files_epoched{1}=[workingdir '/espm8_meg1.mat'];

if(1),
    if(0)
        datadir=[tilde '/homedir/vols_data/murphy_project/data']; % this is the directory the oat files will be stored in
    else
        datadir=[testdir '/faces_group_data_new_full'];
    end;

    spmfilesdir=[datadir '/spm_files']; % this is the directory the SPM files will be stored in

    spm_files_epoched{1}=[spmfilesdir '/efspm8_nosss_meg15_1.mat'];
    D=spm_eeg_load(spm_files_epoched{1});
end;

%%%%%%%%%%%%%%%%%%%
%% DO REGISTRATION AND RUN FORWARD MODEL BASED ON STRUCTURAL SCANS
% This can just be done inside the oat.source_recon stage
% However, it is worth running separately to check that
% the results look reasonable!

for i=1:length(spm_files_epoched),

    %D=spm_eeg_load(spm_files_epoched{i});
    %fidnew=D.fiducials;
    %fidnew.fid.label

    S2=[];

    S2.fid_label.nasion='Nasion';        
    S2.fid_label.lpa='LPA';
    S2.fid_label.rpa='RPA';

    %S2.D = spm_files{i};    % requires .mat extension    
    S2.D = spm_files_epoched{i};    % requires .mat extension    
    
    S2.mri=structural_files{i}; % set S2.sMRI=''; if there is no structural available        
    S2.useheadshape=1;

    %S2.sMRI=S2.mri;

    S2.forward_meg='Single Shell';
    %S2.forward_meg='MEG Local Spheres';

    S2.neuromag_planar_baseline_correction='vector_view';

    D=osl_forward_model(S2);

end;

ca;
%spm_eeg_inv_checkmeshes(D);
spm_eeg_inv_checkdatareg(D);
%spm_eeg_inv_checkforward(D, 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SETUP BEAMFORMER AND FIRST-LEVEL GLM OAT
% In this section we will do a wholebrain beamformer, followed by a trial-wise
% GLM that will correspond to a comparison of the ERFs for the different
% conditions.

%% SETUP THE OAT:
oat=[];
%oat.source_recon.D_continuous=spm_files;
oat.source_recon.D_epoched=spm_files_epoched;
oat.source_recon.conditions={'Motorbike','Neutral face','Happy face','Fearful face'};
%oat.source_recon.freq_range=[0 48]; % frequency range in Hz
%oat.source_recon.time_range=[-0.2 0.4];
oat.source_recon.time_range=[-0.1 0.3];
oat.source_recon.method='beamform';
%oat.source_recon.method='mne';
%oat.source_recon.method='champagne';

oat.source_recon.gridstep=12; % in mm, using a lower resolution here than you would normally, for computational speed
oat.source_recon.mri=structural_files;
oat.source_recon.dirname=[spm_files{1} '_wideband_beamform_sss_PLANMAG'];
%oat.source_recon.modalities={'MEGPLANAR','MEGMAG'};
oat.source_recon.modalities={'MEGPLANAR'};
oat.source_recon.pca_dim=-1;oat.source_recon.force_pca_dim=0;
%oat.source_recon.pca_dim=50;oat.source_recon.force_pca_dim=1;

oat.source_recon.forward_meg='MEG Local Spheres';
%oat.source_recon.forward_meg='Single Shell';

oat.source_recon.work_in_pca_subspace=1;

% design_matrix_summary is a parsimonious description of the design matrix.
% It contains values design_matrix_summary{reg,cond}, where reg is a regressor no. and cond
% is a condition no. This will be used (by expanding the conditions over
% trials) to create the (num_regressors x num_trials) design matrix:
design_matrix_summary={};
design_matrix_summary{1}=[1 0 0 0];design_matrix_summary{2}=[0 1 0 0];design_matrix_summary{3}=[0 0 1 0];design_matrix_summary{4}=[0 0 0 1];
oat.first_level.design_matrix_summary=design_matrix_summary;

% Alternatively you can specify trial-wise subject-specific design matrices
% as ASCII files
if(0),
    % Need to create design matrix without rejects. D.pickconditions auto
    % excludes rejected trials:
    De=spm_eeg_load(spm_files_epoched{1});
    Ntrials=De.ntrials;
    X=zeros(Ntrials,length(oat.source_recon.conditions));
    tot=0;
    for ii=1:length(oat.source_recon.conditions),
        num=length(De.pickconditions(oat.source_recon.conditions{ii}));
        X(tot+1:tot+num,ii)=1;
        tot=tot+num;
    end;
    X=X(1:tot,:);
    
    % save as ASCII file:
    save([workingdir '/subject1_desmat.txt'],'X','-ascii');
    design_matrix_summary={};
    design_matrix_summary{1}=[workingdir '/subject1_desmat.txt'];
    
    try,rmfield(oat.first_level,'trial_rejects');catch, end;
    oat.first_level.design_matrix_summary=design_matrix_summary;
end;

do_hmm=0;
if(do_hmm)
    oat.source_recon.hmm_num_states=-1;
    oat.source_recon.hmm_num_starts=2;
    oat.source_recon.hmm_pca_dim=30;
end;

% contrasts to be calculated:
oat.first_level.contrast={};
oat.first_level.contrast{1}=[3 0 0 0]'; % motorbikes
oat.first_level.contrast{2}=[0 1 1 1]'; % faces
oat.first_level.contrast{3}=[-3 1 1 1]'; % faces-motorbikes
oat.first_level.contrast{4}=[0 -1 0 1]'; 
oat.first_level.contrast_name={};
oat.first_level.contrast_name{1}='motorbikes';
oat.first_level.contrast_name{2}='faces';
oat.first_level.contrast_name{3}='faces-motorbikes';
oat.first_level.contrast_name{3}='fearful-neutral';
oat.first_level.cope_type='acope';
oat = osl_check_oat(oat);

%% RUN THE OAT:
oat.to_do=[1 1 0 0];
oat.do_plots=1;
oat.first_level.name=['wholebrain_first_level'];
oat = osl_run_oat(oat);

res=osl_load_oat_results(oat,oat.source_recon.results_fnames{1});

% load GLM result
stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

% view the GLM design matrix (NOTE that column 1 is motorbikes, columns 2-4 are faces)
figure;imagesc(stats.x);title('GLM design matrix');xlabel('regressor no.');ylabel('trial no.');

% results = osl_get_recon_timecourses( osl_load_oat_results(oat,oat.source_recon.results_fnames{1});, mask_fname )
% osl_oat_plot_hmm_states( oat )

%% OUTPUT SUBJECT'S NIFTII FILES
% Having run the GLM on our source space data, we would like to inspect the
% results for our single subject. 
% We can do this by saving the contrast of parameter estimates (COPEs) and 
% t-statistics for each of our contrasts to NIFTI images.

S2=[];
S2.time.reduce_time=1;
S2.resample_method='fft_ds'
S2.time_ds=2;
S2.time_range=[-0.1 0.3];
S2.data=osl_load_oat_results(oat,oat.first_level.results_fnames{1});
[stats] = osl_reduce_data_to_visualize(S2);

S2=[];
S2.oat=oat;
S2.stats=stats;
S2.first_level_contrasts=[3]; % list of first level contrasts to output
S2.resamp_gridstep=oat.source_recon.gridstep;
[statsdir,times,count]=osl_save_nii_stats(S2);

% VIEW NIFTII RESULTS IN FSLVIEW
% We can now view the nifti images containing our GLM results in FSL, here
% we are runnin fslview from the matlab command line, but you do not need
% to - you can run it from the UNIX command line instead.

mni_brain=[OSLDIR '/std_masks/MNI152_T1_' num2str(S2.resamp_gridstep) 'mm_brain']; 

% INSPECT THE RESULTS OF A CONTRAST IN FSLVIEW. Recall:
%S2.contrast{1}=[3 0 0 0]'; % motorbikes
%S2.contrast{2}=[0 1 1 1]'; % faces
%S2.contrast{3}=[-3 1 1 1]'; % faces-motorbikes

con=3;
runcmd(['fslview ' mni_brain ' ' [statsdir '/cope' num2str(con) '_' num2str(S2.resamp_gridstep) 'mm'] ' ' [statsdir '/tstat' num2str(con) '_' num2str(S2.resamp_gridstep) 'mm']  ' ' [statsdir '/tstat' num2str(con) '_mip_' num2str(S2.resamp_gridstep) 'mm'] ' &']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Investigating locations/regions of interest using a whole-brain OAT
% In this section we will interrogate the wholebrain OAT (run above) using 
% MNI coordinates.

oat.source_recon.dirname=[spm_files{1} '_wideband_beamform'];
oat.first_level.name='wholebrain_first_level';
oat=osl_load_oat(oat.source_recon.dirname,oat.first_level.name);

stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});
% Set the mni_coord to a location of interest
mni_coord=[38,-64,-14]; 

% Find the nearest index of the beamformed voxels to the specified mni_coord 
dists=(sqrt(sum((stats.mni_coord-repmat(mni_coord,length(stats.mni_coord),1)).^2,2)));
[dist,vox_coord]=min(dists); % vox_coord is the voxel index of the beamformed voxels

% Load in the stats structure outputted by OAT first level trial-wise GLM analysis.
%stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});
  
% Plot the parameter estimates for our closest beamformed voxel over time, for 
% our two contrasts of interest:
figure;subplot(1,2,1);
con=1; errorbar(stats.times, squeeze((stats.cope(vox_coord,:,con))),squeeze((stats.stdcope(vox_coord,:,con))),'r'); hold on;
con=2; errorbar(stats.times, squeeze((stats.cope(vox_coord,:,con))),squeeze((stats.stdcope(vox_coord,:,con))),'g');
con=3; errorbar(stats.times, squeeze((stats.cope(vox_coord,:,con))),squeeze((stats.stdcope(vox_coord,:,con))),'b');
legend('Motorbikes','Faces','Faces-Motorbikes','Location','NorthWest');
a=axis; axis([-0.2 0.4 a(3) a(4)]);
plot4paper('time (s)',[oat.first_level.cope_type]);
title([oat.first_level.cope_type]);

% Plot the 1-tailed t-stats over time for each contrast:
subplot(1,2,2);
con=1; plot(stats.times, squeeze((stats.cope(vox_coord,:,con))./stats.stdcope(vox_coord,:,con)),'r','LineWidth',2); hold on;
con=2; plot(stats.times, squeeze((stats.cope(vox_coord,:,con))./stats.stdcope(vox_coord,:,con)),'g','LineWidth',2); 
con=3; plot(stats.times, squeeze((stats.cope(vox_coord,:,con))./stats.stdcope(vox_coord,:,con)),'b','LineWidth',2); 
legend('Motorbikes','Faces','Faces-Motorbikes','Location','NorthWest');
a=axis; axis([-0.2 0.4 a(3) a(4)]);
plot4paper('time (s)','1-tailed t-stat');
title([oat.first_level.cope_type]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Investigating locations/regions of interest using a ROI OAT
% This section will re-rerun just the first level ERF analysis using a ROI 
% in the temporal occiptial fusiform cortex

% Load the wholebrain OAT (which was run above), to make use of the settings 
% and source_recon results already in there
oat.source_recon.dirname=[spm_files{1} '_wideband_beamform'];
oat.first_level.name='wholebrain_first_level';
oat=osl_load_oat(oat.source_recon.dirname,oat.first_level.name);

% Give the first level analysis a new name to avoid copying over whole
% brain analysis:
oat.first_level.name='roi_first_level'; 

% Provide the roi mask as a niftii file:
oat.first_level.mask_fname=[OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex'];
 
oat = osl_check_oat(oat);

oat.first_level.cope_type='acope';

oat.to_do=[0 1 0 0];

oat = osl_run_oat(oat);

%% Spatially average the results over the ROI and plot the results

% load GLM result
stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

% Spatially average
S2=[];
S2.oat=oat;
if(strcmp(oat.first_level.cope_type,'cope')),
    S2.do_resolve_sign_ambiguity=1; 
end;
S2.stats_fname=oat.first_level.results_fnames{1};
if(1),
    S2.mask_fname=[OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex'];
else
    % note that you can use a single coordinate (or a list of coords) as
    % well):
    mni_coord=[38,-64,-14]; 
    S2.mask_fname=osl_mnicoords2mnimask(mni_coord,1,'mask');
end;    
[stats,times,mni_coords_used]=osl_output_roi_stats(S2);

% Plot the parameter estimates for our closest beamformed voxel over time, for 
% our contrasts of interest:
figure; subplot(1,2,1);
con=1; errorbar(stats.times, squeeze((stats.cope(1,:,con))),squeeze((stats.stdcope(1,:,con))),'r'); hold on;
con=2; errorbar(stats.times, squeeze((stats.cope(1,:,con))),squeeze((stats.stdcope(1,:,con))),'g');
con=3; errorbar(stats.times, squeeze((stats.cope(1,:,con))),squeeze((stats.stdcope(1,:,con))),'b');
a=axis; axis([-0.2 0.4 a(3) a(4)]);
plot4paper('time (s)',[oat.first_level.cope_type]);
title([oat.first_level.cope_type]);
legend('Motorbikes','Faces','Faces-Motorbikes','Location','NorthWest');

% Plot the 1-tailed t-stats over time for each contrast:
% Note that these t-stats will be a bit inflated due to the pooling over
% voxels not taking into account the spatial correlation. However, strict 
% statistical validity is not a major concern at the first level; we worry
% about this instead at the group level (and use permutations stats)
subplot(1,2,2);
con=1; plot(stats.times, squeeze((stats.cope(1,:,con))./stats.stdcope(1,:,con)),'r','LineWidth',2); hold on;
con=2; plot(stats.times, squeeze((stats.cope(1,:,con))./stats.stdcope(1,:,con)),'g','LineWidth',2); 
con=3; plot(stats.times, squeeze((stats.cope(1,:,con))./stats.stdcope(1,:,con)),'b','LineWidth',2); 
plot4paper('time (s)','1-tailed t-stat');
a=axis; axis([-0.2 0.4 a(3) a(4)]);
title([oat.first_level.cope_type]);
legend('Motorbikes','Faces','Faces-Motorbikes','Location','NorthWest');

% Now try oat.first_level.cope_type='cope'; 
% Note that we can use the oat.to_do=[0 1 0] to avoid unnecessarily
% rerunning the beamformer

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ROI Time-frequency analysis
% This section will re-rerun the first level analysis using a ROI in 
% the temporal occiptial fusiform cortex, and do a time-frequency
% trial-wise GLM first-level analysis

% Load the wholebrain OAT (which was run above), to make use of the settings 
% and source_recon results already in there
oat.source_recon.dirname=[spm_files{1} '_wideband_beamform'];
oat.first_level.name='wholebrain_first_level';
oat=osl_load_oat(oat.source_recon.dirname,oat.first_level.name);

% Give the first level analysis a new name to avoid copying over previous
% first-level analyses:
oat.first_level.name='roi_tf_first_level'; 

oat.first_level.tf_freq_range=[1 40]; % frequency range in Hz
oat.first_level.time_range=[-0.2 0.4];
oat.first_level.tf_num_freqs=12;
oat.first_level.tf_method='hilbert';
oat.first_level.space_average=0;
oat.first_level.tf_hilbert_freq_res=8;
oat.first_level.mask_fname=[OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex'];

oat = osl_check_oat(oat);

oat.to_do=[0 1 0 0];
oat = osl_run_oat(oat);

%% Spatially average the results over the ROI and plot the results

% load GLM result
stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

% Spatially average
S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.mask_fname=[OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex'];
   
[stats,times,mni_coords_used]=osl_output_roi_stats(S2);

% Plot the parameter estimates for our closest beamformed voxel over time, for 
% our two contrasts of interest:
figure;
con=3; imagesc(stats.times, stats.frequencies, squeeze((stats.cope(1,:,con,:)))');axis xy;
ylabel('frequency (Hz)'); xlabel('time (s)'); colorbar; title(['cope' num2str(con)]);
figure;
imagesc(stats.times, stats.frequencies, squeeze(stats.cope(1,:,con,:)./stats.stdcope(1,:,con,:))',[0 20]);axis xy;
ylabel('frequency (Hz)'); xlabel('time (s)'); colorbar; title(['tstat' num2str(con)]);

% Plot the time course of a single freq bin:
figure;
freqbin=nearest(stats.frequencies,8);
con=1; plot(stats.times, squeeze((stats.cope(1,:,con,freqbin))./stats.stdcope(1,:,con,freqbin)),'r','LineWidth',2); hold on;
con=2; plot(stats.times, squeeze((stats.cope(1,:,con,freqbin))./stats.stdcope(1,:,con,freqbin)),'b','LineWidth',2); 
con=3; plot(stats.times, squeeze((stats.cope(1,:,con,freqbin))./stats.stdcope(1,:,con,freqbin)),'g','LineWidth',2); 
legend('Motorbikes','Faces','Faces-Motorbikes');
xlabel('time (s)'); ylabel('1-tailed t-stat'); title([num2str(stats.frequencies(freqbin)) ' Hz']);

% try comparing hilbert to morlet
% Note that we can use the oat.to_do=[0 1 0] to avoid unnecessarily
% rerunning the beamformer
